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Abstract 


Nontrivial benchmark solutions are developed for the galactic heavy- 
ion transport equations in the straight-ahead approximation with energy 
and spatial coupling. Analytical representations of the ion fluxes are 
obtained for a variety of sources with the assumption that the nuclear 
interaction parameters are energy independent. The method utilizes 
an analytical Laplace transform inversion to yield a closed-form repre- 
sentation that is computationally efficient. The flux profiles are then 
used to predict ion dose profiles, which are important for shield-design 
studies. 


1. Introduction 

With our ever- increasing interest in establishing 
mankind’s presence in space, the protection of per- 
sonnel from energetic radiation in the space envi- 
ronment has become a relevant concern. As high- 
energy radiation interacts with target nuclei, the ions 
undergo nuclear fragmentation and energy degrada- 
tion. This radiation is composed of heavy ions, called 
galactic cosmic rays (GCR), originating in deep space 
and/or protons from the Sun. The nuclear frag- 
ments generated by direct nuclear impact or electro- 
magnetic dissociation form a secondary radiation 
field that again interacts with the target nuclei and 
fragment to add to the biological radiation hazard 
represented by the incident ions. Thus, to ensure 
that the habitat in the space environment is properly 
shielded from energetic radiation, a knowledge of the 
changing nature of the incident radiation field as it 
penetrates protective spacecraft shielding is imper- 
ative. To anticipate future space shielding require- 
ments, NASA has initiated an effort to formulate 
computational methods that simulate radiation ef- 
fects in space. Deterministic transport algorithms 
have been developed for the Boltzmann equation that 
describe GCR interactions with nuclei and cause the 
inevitable changing composition of the incident ra- 
diation field while traversing bulk shielding mate- 
rial. An important component of the NASA pro- 
gram is the assessment of the accuracy of proposed 
deterministic algorithms. For this reason, analyti- 
cal benchmark solutions to mathematically tractable, 
galactic cosmic ray problems have recently been de- 
veloped (refs. 1 and 2). Even though these prob- 
lems involve some simplifying assumptions of the as- 
sociated GCR physics, they still contain the essential 
features of the transport processes. These solutions, 
when compared with the corresponding numerical re- 
sults from algorithms, provide assurance of proper 
programming and a measure of the accuracy of the 
numerical methods used in the algorithm. The moti- 


vation for such comparison is that algorithms devel- 
oped for realistic problems must also yield reliable 
results for the less complicated problems. 

In this report, the analytical solution to the galac- 
tic ion transport (GIT) equations that describe the 
straight-ahead motion of energetic heavy ions (HZE) 
is developed for the first time. An analytical rep- 
resentation of the ion fluxes for a variety of sources 
is obtained with the assumption that the nuclear in- 
teraction parameters are energy independent. The 
method utilizes an analytical Laplace transform in- 
version which yields a closed-form representation 
that is also computationally efficient. From the flux 
profiles, it is then possible to obtain ion dose profiles, 
which are important for shield design. 

In section 2, the analytical solution technique is 
presented in detail, along with solutions for several 
source distributions. In section 3, the numerical 
implementation of the flux and dose profiles is dis- 
cussed. Results for selected case studies and com- 
parisons are presented in section 4, with special em- 
phasis on understanding the behavior of the ion flux 
from both the physical and mathematical points of 
view. The report concludes with a summary and rec- 
ommendations for future development. 

2. Theory 

2.1. Galactic Ion Transport Equations 

Because of the high energy of the galactic ions 
(refs. 3 to 5), the straight-ahead approximation can 
be introduced into the Boltzmann equation with a 
high degree of confidence. In this approximation, 
ions are assumed to move without angular deflection; 
and, as the colliding ions break up in nuclear frag- 
mentation, the fragments continue in the incident ion 
direction. Thus, for ions of charge number j, the 



appropriate transport equation is (ref. 6) 


d_ 

dx 




<pj{x,E) 


J 

= M lk {E)<j k {E)4> k {x,E) (la) 
k=j + 1 


where <f>j is the flux of the jth ion flux at position x 
with an energy per nucleon of E. The macroscopic 
absorption cross section is , Mj is the multiplicity 
of ion j produced in collision with ion k , and Sj is 
the absolute value of the change in energy E per 
unit distance traveled. The added assumptions of 
continuous slowing-down theory and that target ion 
fragmentation can be neglected have been made in 
equation (la). Source ions are assumed incident on a 
semi-infinite shield with an initial energy distribution 
of 

4 j (0,E) = fj{E) (lb) 

Because of the success of the analytical investigations 
applied to previous benchmarks (refs. 1 and 2), the 
analytical approach is again followed here. To make 
equation (la) mathematically tractable, the cross 
sections and multiplicities are assumed to be energy 
independent. This is a reasonable assumption at high 
energy. Also, the well-known scaling law for stopping 


powers at high energies (ref. 6), or 


Sj(E) = Vj S P (E) 

(2a) 

is introduced into equation (la), where 


vj = Zj /Aj 

(2b) 


and where S p is the proton stopping power and Aj 
is the atomic mass of the jth ion. With this substi- 
tution, and with the change of variable from energy 
to proton path length at a given energy s(E, Eq), 
equations (1) become 


d d 
]h: +Uj ds +aj 


4>j{x,s)= ^2 0j k <t>k(x,s) 
k=j+\ 


and 

where 


4>j{0,s) = fj(s ) 


Pjk - Mjk a k 

The path length is given by 

r^o 


» (E, Eq) = [ ° dE'/S p (E') 
J E 


(3a) 

(3b) 

(3c) 

(3d) 


which, from the definition of the flux distribution, 
implies 

4>j(x,s) = Sp(E) 4>j(x,E) (3e) 

and 

f j (s) = S p (E)f j (E) (3f) 


Some readers may be accustomed to seeing the en- 
ergy variable in equations (la) and (lb) transformed 
by using the residual range R p {E) (refs. 7 and 8), 
where 



dE f 
S{E> ) 


The variables R p {E) and s(E,Eq) are directly 
related. The variable s{E, Eq) is the penetration 
distance for a proton of initial energy Eo to reach 
present energy E\ R{E) is the remaining distance 
the proton will travel before stopping. 


The major difficulty in applying the path-length 
transformation is in the determination of s(E , Sq) 
from the energy E by using equation (3d), since an 
integral must be inverted to return to the energy- 
dependent flux given by equation (3e). This added 
inconvenience, however, is readily acceptable be- 
cause of the enormous simplification offered by the 
energy path- length transformation. Indeed, equa- 
tions (3a) and (3b) now closely resemble the energy- 
independent case considered in reference 1. This sim- 
ilarity suggests that similar solution methods should 
therefore be adopted. 


While the emphasis here is on generating highly 
accurate benchmark solutions, the only limiting as- 
sumption made to obtain equations (3) is that of 
constant cross sections and multiplicities. Thus, the 
analytical results to be presented can be of use as 
a predictive as well as a benchmarking tool to give 
insight into realistic shield design. 


2.2. Analytical Solution for 
Monoenergetic Source 

Initially, the solution will be obtained for a mono- 
energetic beam source of ions of charge J and en- 
ergy Eq(s = 0) incident on the semi-infinite shield 
lfj{E) — S(E — Eq) 6ij]. For this case, <f>j corre- 
sponds mathematically to the Green’s function, since 
the source is of the form 


<t>j{0,s) = 6{s) 6ij (3g) 

which translates into the following volume source: 
,Q{x,s) = 6(s) 6(x) 8 U 

The Green’s function will be used subsequently to 
generate solutions for more general sources, includ- 
ing beams composed of several different ion types 


2 


(composite beams). For the incident ion J and for 
ion J - 1, produced by target fragmentation, equa- 
tions (3) can be solved analytically. It can be shown 
that the solutions for ions J and J - 1 are as follows: 


where the transform is defined by 



ds exp(— sp) (p(x, 


«) 


< t>j(x , s) = exp(-ajx) 6{s - vjx) (4a) 


= ^-^{exp (-ojx) 

x exp[-ai(s - vjx)/b\^Q{s - vjx) - exp(-ffj_ix) 
x expf-ai(s - vj-\x)/bA 0(s - 

^ (4b) 

where a\ = h\ = — and 6(a) is the 

Heaviside step function. The increased complexity of 
these solutions relative to the benchmarks considered 
in reference 2 is clearly evident. Because of the 
interaction of x and s through the transport operator, 
a wave of particles in x, s space is induced. The 
delta function in the expression for <j>j represents that 
wave for the incident ions. The delta function is a 
mathematical statement required by the continuous 
slowing-down approximation. For ions of type J to 
be at a given position x, measured from the surface, 
they must have traversed a path length sjvj exactly. 
The leading exponential represents the attenuation 
(loss) of these ions as they fragment into lighter ions 
along their trajectories. Ion J — 1, produced by the 
fragmentation of ion J, cannot exist before its source 
(at x = s/vj). Also, because of the finite range of 
a charged particle, ion J — 1 cannot travel farther 
than x = s / v j _\. Thus, <£j_ \ is nonzero only in 
the interval sjvj < x < sjvj_\, as indicated by 
equation (4b). This technique is better illustrated 
in the (x,s) space diagram in figure 1, where the 
region of nonzero solution is indicated for each ion. 
The above analysis holds for all ions j produced 
by fragmentation, and the farthest nonzero solution 
boundary approaches s = x as j goes to one. The 
region x > s is forbidden, since no ion has yet 
traversed a path long enough to get there. 

At first glance, one would not attempt to carry on 
the analytical manipulation necessary to determine 
the flux for the lower values of j. Fortunately, the 
analytical results already derived in reference 2 are 
very useful in pursuing the solution here. 

Taking a Laplace transform of s in equations (3a) 
and (3b) gives 



0j(O,p) = &ij (5b) 


Treating equations (5a) and (5b) as a set of ordinary 
differential equations yields the general solution 

i 

4>j_ i (x,p) - a i,J-e(p) cx Pl^( (T J-i + v J-iP) x \ 
1=0 , v 

(6a) 

where the coefficients depend on the trans- 

formed variable p and are given recursively by (ref. 2) 


l 

°i,j-r(p) - _ <Tj _ i ) - p(vj_ ( - vj_i ) 

(-i 

x 0J-U-V OiS-eV) (0<^^D 

k'—i 

(6b) 

and 

t-1 

<*e,j-dp ) = a i,J-e(P) ( 6c ) 

t=0 

The last expression is necessary to satisfy the bound- 
ary condition <fij(Q,p) = 0 for 1 < j < J — 1* 

At this point, a numerical inversion procedure 
could be attempted. Such an attempt, however, 
would most likely fail, since it is known a priori that 
the expression for contains step functions that 

characterize the spatially discontinuous ion density 
waves; and in general, numerical inversions fail near 
such discontinuities. Thus, a numerical inversion will 
not work well in this case, even for the heavy ions. 
Therefore, the only other option is to perform the 
inversion analytically, as was done for ions J and J— 1 
(cqs. (4a) and (4b)). To do so, however, requires the 
explicit knowledge of the dependence of c*i y j-£ on p, 
which will now be determined. 

For ion J — 1, we have 


«o.j(p) = 1 


a 0,J-l(p) = 


aij-i(p) = 


r 0,l 

+ P b J-i,J 

\r 1 

*1,1 

a J-i,J +pbj-i,j J 


(7a) 
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where 


Ur' = * u ' 

Yu'=0J- l,J 

> 

a £m — a l ~ &7n 


« 2 „/- 2 (p) = v 2 y 2 A,(p) + y 2 J 2 2 A 2 fp) + Y 2 J - 2 A 3 (p) (10b) 

where 


(7b) 


A 3 (p) = 


1 

a J—2,J—l + P b J-2,J-\ 


^hri — 17 £ 17 m J 


For £ — 2 and i — 0, 


It should be noted that Ax is common to both the 
£ = 1 and t — 2 a-coefficients. A condensed expres- 
sion for oli j _ 2 is therefore 


«<),./- 2 (?) = 


a J-2,J + Pbj-2,j 


\ a J~lJ +pbj-\ i j 
Upon using the partial fraction expansion 

1 1 

a J-2,J +pbj-2,j aj-\,j +pbj-i t j 

_ ^J-2,J 

a J-2,J b J-l,J ~ aj-i,j bj_ 2 ,j 

+ 6 J~l ,J 

a J-2,J bj- 1,J ~ a J-l,J bj-2,j 


+ 0J-2,J (8) 


Ai(p) 

A 2 (P) 


3 

«U-2(p) = E^r" 2 Ar(p) (11) 

r — 1 

with the Y coefficients set to zero as appropriate. 
When the last equation is introduced into equa- 
tions (6b) and (6c) with j = J - 3, the following 
form for ot% is found: 

6 

“i,J-3(p) = E Y i~ 3A r(p) (0 <i<?) 

T — 1 

where the rule for expansion in partial fractions 


where 


and 


AlC P) = 


1 

+ pbj-ij 


A 2 (p) = 


1 

a J—2,J + Pbj-2J 


we obtain 


£^m{p)A-n{p) ^1 mn^m(p) "I" ^2 mn^n (p) 

has been used. Apparently, with each unit increase in 
£ ) £, new terms appear in the expression for a ± j_^(p), 
with the appropriate Y coefficients again set to zero. 
Thus, the general expression for o^j_^(p) is 


«0,J-2(p) = Vfff 2 Ai(p) + V'^- 2 A 2 (p) 
with 

yj - 2 = Pj-\J 

1)1 a J-2>J b J- IJ - a J-l,J b J- 2J 

and 


(9a) 


(9b) 


<*i,J-e(P ) = E n>"' A riP ) M = t(£ + l)/ 2 ) 

r= 1 

( 12 ) 

and the explicit dependence of a j j-t(p) on p has 
been established. Table I shows how, given an 
index r, the appropriate values of a m ^ n and 6 m>n for 
the term 


\rJ—2 

1,2 ' 


Pj- 


2 J 


b J-2 t J Pj-2,J-l 
°J— 2, J b J-\,J ~ a J-\ } J b J-2J 


(9c) 


The coefficients of A r (p) are independent of p. Simi- 
larly, for a h j_ 2 {p) and a 2 j- 2 (p), 


oi,7-2(p) = >3f 2 Ai(p) + n3“ 2 A3(p) 


(10a) 


and 


A rip) 


1 

a m,n T bm t nP 


are determined. A general recurrence relation for 
Y i r 1 is obtained by introducing equation (12) into 

the relation for aj j_^(p) (eqs. (6b) and (6c)) to yield 
(after much algebra) 
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J—i 

k = J-t+ 1 

(l<r<JV*_ lt 0<i<e-l) (13a) 

J-i N J-k 

Y i/ = E £/-/,* Y, *ij!r Y £,r + 0J-t,jKo 
k=J-i+\ r=l 


The inversion is easily performed analytically, since 


r-l 


ex P(-^y 


-i*P)l 1 r Or. J 

~J = 6 ; exp [-*<— *-**>1 


9(s - I'y-ji) 


Or + b r p | 6 r " *' [ b r 

to give the following solution for 1 < i < J - 1: 

' N t Y J-l 

<t>J-e{x , «) = 51 exp(-<r i/ _ i a:)0(s - vj_ iX ) 


i=0 


(r* = N(_i + i + 1, 0<i<l-l) (13b) 


Y. 

IT 


J-i 


= 0 


( N t-l + 1 <r< N e _ x + i, N(_ l +i + 2<r<N e ) (13c) 


e - 1 




i=0 


(1 < r < Ay i = t) (13d) 


where 


= _ 

Mi. r 


vj-itJ-i 


-a r bj_g j_ i + j b r 


X6XP (15) 

Thus, a closed-form analytical solution to the GIT 
equation for constant properties has been derived. 
A formulation more fundamental and therefore more 
appropriate for computational purposes is obtained 
from equation (15). 

As mentioned previously, some of the coefficients, 
Y i r , are identically zero in order to write a z j-t{p) 
as the simple summation of equation (12). This pro- 
cedure dramatically increases the required computer 
storage and also increases the potential for the ac- 
cumulation of roundoff error. These zero coefficients 
can be conveniently eliminated by noting the physical 
requirement that 


*2,i , r — 


b r 


a r — a m ^ n ^ 


<f>j-i(x,s) = 0 (16) 

for x < s/i/j and by interchanging the summations 
in equation (15) to give 


br — ^m(r),n(r) 


and 


h(m) = 


f 1 m > 0 
l 0 m < 0 


4>j-((x y s) 


N ( 

E 


exp(-^s) 


br 


E 


x exp - ^i//_,-)x| Y t J r 1 8{s - vj-ix) (17) 


Interchanging the summations in equation (13b) 
yields the following further simplification: 


N t J-i 

Y i,r* ~ 5 . x l.«.r ^ Pj-e,k Y i t r b{Nj_ k - r) + 0 J-(,J Sift 
r=l *=•/-*+! 


N i-i /h n 

■sw-f 


+ Pj-i,j Ao 


(13e) 


The final expression for the transform is obtained 
by substituting equation (12) into equation (6a) as 
follows: 


t N e 

4>J-i( x ,P) = ^2 exp(-<y t/ _jg) 

i-0 r~l 


Y: 


J-i 


ex p(-^-t^p) 
a r + b r p 


(14) 


Equation (16) comes about because ion fragments 
cannot be produced before the incident ion has trav- 
eled to its position (x = sjvj ) as mandated by the 
continuous slowing-down theory. Since in most cases 
of interest, vj > Vj _ x for 1 < i < l, we have 
4>{x,s) = 0, where s > xvj > xvj_ { . Also from 
equation (17), 


N { 

° = E 


r=l 



£ 


E 


x exp 


K 


^ ~ V r Vj ~ l 1 * 


y 


j-t 


l,T 


(18) 


which must hold for all values of x and s that 
satisfy the above inequality. Equation (18) is trivially 
satisfied for x = 0 as a result of equation (13d) of the 
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recurrence relations for Y^ T ^ ■ For equation (18) to 
be satisfied in general, the following must hold at 
each value of r and i, independently of x : 



CLf 

W 


VJ-i)x 


Y i J r~ e = 0 
i.r 


(19) 


One way for equation (19) to be a true statement is 
for the coefficient of x to be identically zero. This 
requirement, however, is obviously too restrictive. A 
second possibility is that the coefficients vanish at 
several values of i taken in pairs; when this is not the 
case, Y/- -£ should be zero. It can be observed that 


aj. H - fvj-i, = cj. h - (20) 


for two values of 1 ( 21 , 12 ) each value of v. Also, as 
required by equation (19), 


yJ—t _ 

H,r hS 


(21) 


advantage of the analytical solution (eq. (22a)) in 
the path-length variable s is that it is independent of 
the energy range relation. Once the proton stopping 
power is specified, then <^j_^(x, E) is obtained from 
equation (3e). Additional quantities of interest, as 
well as solutions to more comprehensive benchmarks 
to be considered in the following sections, arc ob- 
tained from the analytical form of Green’s function. 

2.3. Dose Profile 

The dose profile (in x) for the ion J — £ (ref. 2) 
is given by integrating the total track length of 
ion j over all energy in a differential volume element 
multiplied by the energy loss per track length as 
follows: 

ro o 

D j-t(x) = / dE Sj—£ (f>j-e(x, E) (24) 

Jo 

Recall that Aj^e is the atomic mass of the J — l 
ion. Upon substituting the stopping-power scaling 
law (eq. (2a)) and equation (3e) into equation (24), 


values of i\ and i 2 for 1 < r < 15 are given in table II, 
from which a recognizable pattern emerges. 

By introducing equations (20) and (21) into equa- 
tion (17), the expression for , s) then becomes 


N f 

A/-<(z- s ) = X/xp(-^ S ) pxp 

r=l 




D J-d x ) = A J-t v J-t 



ds S p (E(s)) <p{x,s) 
(25) 


where R Pu is the proton range in the shield material 
for incident ions of energy E 0 . W1 icn equation (22a) 
is introduced into the integral, the following closed- 
form representation for the dose results: 


x - » v-t,*) - 0 ( s - 'v-.j*)] 

(22a) 

Dj- 

where 



Y J ~ ( = Y J ~Jlb T 

il,r i\ X f 1 

(22b) 

where 

The recurrence relation for becomes 






K/ = xti E y? u r h (Nj_ k ■ 

-r) 


+ 

1 

II 

(23a) 

Ir (x) 

for 1 < r < AT^_j and 



j 

k,o 


(23b) 



N ( 


r= 1 


£r — <y j~i } — (3 r 
fir — d r /b r 


(26b) 


0 (E) cxp(-j3 r s) 


fv J- ti* 


■0{r o -uj_ i2 x) / ds S p (E)exp(-(3 r s) (26c) 
Jvj-i 2 x 

r 0 = mm(Rp 0 ,isjx) (26d) 


for + 1 < r < N£, where the sum in equa- 

tion 23(b) is over the r' , for which Y % J ~J gives a 
nonzero contribution. A sample of these particular 
values of r' and the parities P r > are given in table III. 

Equation (22a) is the expression for the analytical 
solution of the GIT equations to be evaluated. One 


To perform the integration either numerically or ana- 
lytically, we are now faced with the task of specifying 
the range energy relation for protons. 

Initially, a simplified range energy relation is spec- 
ified so that the integrals I r can be performed analyt- 
ically to guarantee reliable benchmark results. This 
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relation takes the general form (ref. 9) 

s ~ R Po~ a o ln(l + a\ E n °) (27) 

The simplest case is for n Q taken as unity, which 
yields a simplified stopping power 

S p = aexp[(s — R Po )0] (28) 


When the stopping power, as represented by 
equation (28), is introduced into the integrand of I r 
(eq. (26c)), the integrals can be performed analyti- 
cally to give 


Ir(x) = 


a exp (R Po P) 
Pr+P 


0(r o -vj-i jX) 


where 


1 

a = 

OL\a 0 



Oto 


The following prescription is used to determine the 
parameters a and (3: 


| exp - 


(p r + p)vj_ h x 


) jexp 


— exp 


-( Pr + P)r a 


e{r 0 - vj-i 2 x)< exp -(P r + P)vj_ i2 x 


exp 


~{Pr + P)Tq 


]} 


(30) 


(a) The more accurate proton stopping-power (called 
the Wilson stopping power) relation (— dE/ds) 
given by equation (27) with (ref. 9) 


ra 0 =1.79 
a G — 500 

«i = 3.66 x 1(T 6 , 


(29a) 


For the more comprehensive stopping-power formula 
with a Q) c*i, and n Q specified by equation (29a), 
the integrals must be performed numerically. This 
is most conveniently done by transforming to the 
energy variable. Then the integrals in equation (26d) 
are of the form 

fE(vj_ t x) 

I r ,i = / dE exp [~/3 r s(E)] (31a) 

JE(r 0 ) 


is used to specify the proton range as follows: 


where s{E) is given by equation (27) and 


R Po = a 0 ln(l +<*!££■) (29b) 

(b) Since, in the energ^ r variable, the simplified stop- 
ping power is 

S p (E) = a + j3E (29c) 

/? can be extracted as the slope 

(3 = [SpiE^ - S P (E 2 )\/(E X - E 2 ) (29d) 


given two values of stopping power S p (Ei) and 
S p (E 2 ) at energies E\ and E 2> respectively. These 
stopping powers are obtained from the more 
comprehensive stopping power referred to in (a) 
above. 

(c) If E = 0 when s = R Po and E — E 0 when s = 0, 
then 

— — ^ (29e) 

0 exp (R Po {3) — 1 ^ 

gives a consistent value of a with (3 and the proton 
range. 


E(s) 


exp [(R Pu 


- s)/a 0 ] - 1 
«i 


\jn 0 


(31b) 


Finally, / r ? ; is evaluated by using a reliable quadra- 
ture scheme. The dose for the uncollided ions, 
£ = 0, is obtained by substituting equation (4a) into 
equation (24) 


f R Po 

Dj(x) = Ajvj / ds S p (s ) <pj{x,s) 

Jo 


= Ajvj exp (-< tjx ) S p {vjx) 8{R Po - vjx) 

(32) 

for any proton stopping power S p . 

2.4. Source With an Energy Distribution 

Since the solution for a beam source is actually 
Green’s function in the energy variable, a source for 
ions of charge J impinging on the shield surface of 
the form 


4> j (0,E) = 6 iJ f(E)6(E 0 -E) (33) 

can be accommodated. An energy cutoff E 0 has 
been imposed so that the path length given by equa- 
tion (3d) will again have meaning. 


7 


The jth ion flux for this source is therefore 


E ) = [ E ° dE' f(E') <pj(x, E') (34) 

' J E 

where the explicit dependence of Green’s function 

on the beam energy E r has been included. The 
appropriate change of variable yields 

E) = <fr*(x, s)/S p (E) (35a) 

and 

f(E) = f(s)/Sp(E) (35b) 

where 

4>*(x,s) = / ds' f{s - s') (35c) 

J JO 

Then, by substituting equation (22a) into equa- 
tion (35c), the closed- form representation of the flux 
from a distributed source is 

N e 

4>*j_ ( {x,s) = £ exp(-e r x) Y/^ r e J r (x,s) (36a) 

r— 1 

where 

f r ° 

J r (x,s) = 0(r„ - uj_ h x) / ds' f(s - s') exp (-&«') 

Jxvj-i x 

- S{T 0 - vj-i 2 x) f ds' f(s - s') exp (-&•«') 
Jxvj. i„ 

(36b) 

and 

r Q — min(s, xv j) (36c) 

Finally, for the incident ions (f = 0) from equa- 
tion (4a), 

cp*j(x,s) = exp (-< t/x) f{s - vjx) 9(s - vjx ) (36d) 

To generate benchmark results, we consider, as an 
example, a source of the form 

/(s) = Q 0 exp (— M 0 s) (37a) 

or in the energy variable 

m = m 

yields 


Jr(x,s) = j exp (-Mo*) (0(ro ~ vj- ix *) 

x {exp [(Mo - Pr)ro] - exp [{M 0 - 0r) } 

- 6(r 0 - vj-i 2 x){exp\(M 0 - 0r)r o ] 

- exp [(Mo - /3r)z*//-t 2 ] }) (37c) 

For the source ions, 

<f)*j(x,s) = Qo exp \-{aj - M 0 i//)x] exp(-Af 0 s) 9(s - vjx) 

(37d) 

A more general class of analytical sources can also 
provide benchmark flux profiles. If a source f(s ) pos- 
sesses an explicit Laplace transform f(p ), then since 
the integral in equation (35c) is of the convolution 
type, we find after application of a Laplace trans- 
form that _ 

~4>){x,p) = ~f(p ) 4>j(x,p) 

The transform of (j>j(x,s) is easily obtained from 
equation (22a) as 

Ne 

4j(x,p) = £ exp (-e r x)Y^~ e 

exp [-(fir + p)vj-i l x] exp[(/? r 4- p)vj-j 2 x ] 1 

Pr + p + P j 

(38) 

Then, upon inversion of equation (38), 

n £ _ ' 

<t>j(x,s) = £ exp(- fr i) 
r=l 

x [exp(-/3 r ^j_ il x)L r {s - vj_^x) 0{s - vj-i x x) 
-exp(-/3ri 'j-i 2 x)L r {s - vj-i 2 x) 9(s - vj_ i2 x)} , 

(39a) 

where _ 

< 39b > 

If the inversion C~ l can be performed analyti- 
cally, then accurate benchmarks can be generated. 
In most cases, however, the transform £annot be 
performed analytically. For a general /(p), it is 
possible to perform the inversion numerically by 
using a recently developed numerical-inversion algo- 
rithm (ref. 10). Since numerous numerical inversions 
would have to be performed for a heavy incident ion 
( J > 10), the required computational time would 
be considerable. For lighter incident ions ( J < 10), 
the above method could be applied for a reasonable 
computational effort. Of course, an even more gen- 
eral class of sources could be accommodated by using 


8 


a quadrature scheme to evaluate J r (x, s ) as was done 
for the dose. 

2.5. Composite Beam With Distributed 

Sources 

Another generalization of the beam source is the 
composite beam. In this case, several ions of different 
charge are allowed to impinge on the shield surface 
with the same source incident energy per nucleon of 
E 0 (same velocities). The GIT equations (eqs. (3a 
and 3b)) then become 

j 

Lj4>j{x,s)= ^2 Pjk <t>k( x > s ) (1 <3<J) (40a) 

1 

Since the GIT equations are linear, the composite beam can be decomposed into its individual ion 
components with a source at x = 0. Thus, 

J 

Lj<f>j (x, s)= Pjk <t>k (*> s ) </>j (0, s) = dj Sij 6(s) 

k=j+l 

J-l r , 

Lj~i4>j~ 1 (x, s) = y /3 jk <j> J k ~ l {x,s) <l>j (0,s) = d/-i S(s) 

k=j+l 

> ( 41 ) 

2 

L 2 <pj(x,s) = y 0 jk <t> 2 k (x,s) 
k=j + 1 

L\ <j>) (x, s) = 0 


where use was made of the requirement that 
<pj~ e '(x,s) = 0 

for j > J — £ f . Therefore, because of the uniqueness 
of the GIT equations, a comparison of equations (40) 
and (42) yields 

J-l 

= X <pj- £ \x,s) (43) 

f=0 

Since (x,s) is Green’s function, the sum in 
equation (41) is obtained from equation (22a). Sim- 
ilarly, for a distributed source, the sum in equa- 
tion (43) is obtained from equation (36a). 


where the superscript indicates the incident ion 
charge. When equations (41) are summed, we 
obtain 


j-i j-i J-l 

Lj y: <^ -<, (x,s)= y Pjk y (42a) 

^=0 k = j + 1 ^'=0 

and 

j-i j-i 

1 (o? 5 ) = ^2 dj-? j-? 

£'=0 t '=0 

J 

= y df Sjj, 6(a) (42b) 

/=! 


4>j(0,s) = d 2 Sj 2 6(s) 
</>](0,s) = d x 6j\6(s) 
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2.6. Distinction Between Neutrons and 

Protons 

As currently written, the GIT equations treat all 
the ions with a mass number of one as protons. They 
all have an associated stopping power and range. To 
make the GIT equations more realistic, a neutron 
component has been added. Neutrons are given the 
charge number zero (j = 0), and the equation that 
describes their straight-ahead neutron motion is 

(jL + <T °) ^ X> M 0.k<?k<t>k(x> s ) ( 44 ) 

where v o = 0 and where it is assumed that protons 
cannot act as a source of neutrons. The absorption 
cross section is taken to be the same as for protons, 
and the multiplicities are assumed to be a fraction 
f n of those for protons as follows: 

M 0 ,i = 0 1 

(45) 

M 0ik = f n M lik (k > 2) J 

The proton multiplicities are then modified by a 
factor 1 — fn- 

The formulation also assumes that the neutrons 
are monoenergetic and do not have scattering inter- 
actions. The inclusion of neutrons here is done only 
for completeness; a more substantial model is cur- 
rently under development. 

3. Numerical Implementation 

3.1. Nuclear Fragmentation Model and 

Stopping Power 

To numerically evaluate the ion flux and dose 
for benchmarking purposes, the following simplified 
nuclear model is assumed: 


<?j = ° P 


(46a) 


( F^T 

(k > j ) 


Mjk < 

1 


(46b) 


u 

(k < j ) 



The choice of <jj is based on nuclear liquid drop model 
considerations. The multiplicities are chosen so that 
charge is conserved in each interaction. The cross- 
section normalization is an input parameter that is 
taken to be representative of an air shield as follows: 

a = 0.01247 cm 2 /g 


Unless otherwise stated, the simplified stopping 
power is used in the computations of this report. The 
values of energy for the determination of / 3 from equa- 
tion (29d) were taken as E\ = 10 MeV/nucleon and 
E 2 = 10 4 MeV/nucleon, which provided a represen- 
tation over the energy interval of interest. 

3.2. Round-off Error 

Because the solution for the ion flux that results 
from an incident beam is an analytical representa- 
tion, discretization and truncation errors are not rel- 
evant to the numerical evaluation. However, a very 
real concern is the error that results from round off. 
The central numerical procedure used to evaluate 
the flux is the recurrence relation for the coefficients 
(eqs. (23)), which is rather involved. It is well- 
known that significant round-off error can accumu- 
late when using recurrence relations. For this rea- 
son, the simplifications described in section 2.1 play 
an important role in round-off error mitigation. 

The major cause of round-off error can be traced 
to the r summation in equation (22a). Table IV 
shows each term calculated in Control Data Corpora- 
tion (CDC) double-precision arithmetic (~29 digits) 
and the corresponding partial sums of the summa- 
tion over r in equation (22a) for x = 16.004 g/cm 2 , 
E = 100 MeV/nucleon, and j — 16. The incident 
ion is nickel (J — 28) at E 0 — 5 GeV/nucleon. The 
summation is remarkable, in that the partial sum os- 
cillates between ilO 11 with a final result of ^10~ 3 . 
Thus, a precision of at least 14 digits is required for 
this computation to have one correct digit. The cor- 
responding single-precision calculation fails to give a 
positive flux. The failure of single-precision arith- 
metic is expected because of the significant loss of 
precision. All computations are performed in double- 
precision arithmetic. Even though round-off error 
will cause failure of the r sum for the light ions pro- 
duced from the heavy incident ions for J > 30, there 
is no apparent ill effect for J < 30. If necessary 
and available, quadruple-precision routines could be 
employed. Round-off error is also expected to limit 
the dose and distributed source evaluations to only 
heavy ions when the integrals I r and J r cannot be 
performed analytically. 

4. Results 

* 

The flux and dose expressions determined in sec- 
tion 2 are numerically evaluated for air shields in 
this section by using the nuclear model and stopping 
power discussed in section 3. 
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4.1. Flux Profile From a Monoenergetic 

Beam Source 

According to the continuous slowing-down ap- 
proximation, the incident ion that impacts the shield 
surface with energy E 0 has lost a known amount of 
energy after having traversed a given path length s. 
Thus, at a particular energy £q, a fraction of the in- 
cident ions given by exp(— ajx) have survived to the 
distance xj — s\(E\)/uj i where s\ is determined 
from equation (27). At this point, crj(t)j(x,s) dx in- 
cident ions/sec per cross-sectional area fragment into 
lighter ions (0 < j < J — 1). Because of their finite 
velocities, the secondary ions of charge J — 1 that 
originate at xj can then propagate to a maximum 
range of xj^i = s\(E\)/vj_i at this energy (fig. 2). 
The flux profile decays with x because there is no 
source of secondary ions between the creation and 
allowable range. The fluxes <}>j{x, s) for several ener- 
gies E\ < #2 < • . • < E§ are displayed in figure 2(b). 
In keeping with the range energy relation, the region 
for which E) is nonzero moves farther into 

the shield with decreasing energy (increasing path 
length). Also, the flux decreases with decreasing en- 
ergy, because the penetration distances are greater 
at lower energies; these greater distances allow more 
of the J — 1 ions to fragment into lighter ions. Fig- 
ures 3(a) and 3(b) show the variation of <fij_ i(x, E) 
with E at a fixed positron x and at several positions 
xi < X 2 < . . . < Xfy. The decrease in the flux at lower 
energies is clearly evident. 

In the region s\{E\)/vj < x < si{E\)/vj_i, 
the secondary ions (J — 1) continuously generate 
lighter ions. Figure 4(a) shows both 4>j-\{x,E) 
and <pj_ 2 {x,E) at energy E\. In the region where 
ions J — 2 are produced {s\fvj < x < 
the flux increases as it should. When the source 
abates, the flux decreases until the ions reach their 
allowable range x,/_2 = s l/^j-2’ At xj_j, the slope 
of (f>j -2 changes dramatically (being infinite), which 
indicates the reduction of the source abruptly to zero. 
Figure 4(b) shows the secondary and tertiary ion 
fluxes for several energies ( E\ < E<i < . . . < E§). 

Considering the lighter ions ( J — 3, J — 4, J — 5) 
with representative fluxes shown in figures 5(a) and 
5(b), we observe similar behavior. In the region 
s\ j uj < x < s\/uj- 1, where ions J — 1 and J — 2 
produce ion J - 3, the flux increases as a result of 
the increasing source. The slope again changes when 
the source from becomes zero. Because of the 

production of ion J — 2 in the region < x < 

s\jvj_ 2 , and because of the attendant competition 
with loss from fragmentation, the flux for ion J - 3 
initially increases but then decreases. In the region 
of no J - 3 ion source ( s/vj~2 < x < si/^j-3), the 


flux is monotonically decreasing. The flux profiles 
become increasingly smooth for the lighter ions. This 
increase is a consequence of the smaller fraction 
that the discontinuous source represents of the total 
source at < i < l - 2 for ion J - L In 

addition, at where the only source is from 

the preceding ion, attenuation has reduced the source 
contribution significantly. However, a kink in the 
(pj_ 5 flux curve can still be observed at xj_ 4. 

When neon ions of energy 1 GeV/nucleon are 
incident on the shield surface, the resulting fluxes 
of the ions produced by fragmentation are shown 
in figures 6(a) to 6(c) for E = 995 MeV/nucleon, 
950 MeV/nucleon, and 100 MeV/nucleon. For the 
energy nearest the source energy, E = 995 MeV/ 
nucleon, the fluxes are a series of waves of almost 
uniform density that originate at x « 0.6 and extend 
to the maximum range of the proton (j — 1) at this 
energy. Since the ions have penetrated only a short 
distance into the shield, significant fragmentation has 
not yet occurred. For this reason, the maximum 
flux for each successive ion is less than that of the 
preceding ion until j — 1. Because all ions for 
which j > 1 can produce protons, the proton source 
starts to become significant and results in an increase 
in the flux near the end of the proton trajectory. 
This increase in fragmentation is clearly evident at 
950 MeV/nucleon; these are indicated by significant 
increases in the fluxes of light ions j = 1,2,3. At 
100 MeV/nucleon, the ions have almost penetrated to 
their maximum ranges (309.75/^ ), and the light ions 
now have the largest fluxes. Figure 6(d) displays the 
three energies and clearly shows the changing nature 
of the particle flux with penetration that results from 
fragmentation. Figures 7(a) and 7(b) show the flux 
energy dependence at several positions and show that 
the energy spectrum broadens with penetration. 

The incident ion energy influences the flux profiles 
significantly, as shown in figures 8(a) and 8(b), where 
profiles at E = 100 MeV/nucleon for neon ions of in- 
cident energies 500 MeV / nucleon and 2 GeV / nucleon 
are shown. The heavy-ion flux is still present for an 
incident energy of 500 MeV/nucleon; however, for 
2 GeV/nucleon, so much penetration has occurred 
that the heavy ions have almost all fragmented and 
left a negligible heavy- ion flux. 

The cross-section normalization has a major effect 
on the profiles. Figures 9(a) and 9(b) show the 
flux profiles at 995 MeV/nucleon for neon incident 
at 1 GeV/nucleon with cross-section normalizations 
of 1(T 4 cm 2 /g and 1.0 cm 2 /g. The normalization 
is proportional to the shield density. Thus, for 
the rarefied shield (a = 10 -4 ), the ions behave 
as expected with only a small amount of light ion 
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production because of the relatively few interactions 
that cause fragmentation. However, for the dense 
shield (a = 1.0), a significant flux of light ions has 
already been generated after a penetration of only 
3 g/cm 2 . 

To demonstrate the general applicability of the 
analytical solution, the flux profiles for E — 4.995 
and 2 GeV/nucleon for nickel incident on a shield 
at 5 GeV/nucleon are presented in figure 10. To 
obtain better resolution for the 2-GeV case, (f>j for 
8 < j < 27 and 1 < j < 7 are plotted on separate 
plates. With GDC double precision, incident ions of 
charge up to ^30 can be considered. For J > 30, 
the round-off error becomes a problem and higher 
precision routines must be used at the expense of 
storage and computational time. 

4.2. Dose Profiles From a Monoenergetic 

Beam Source 

The dose profiles given by equations (26a) to 
(26d) for the stopping power of equation (27) are 
evaluated here. Figures 11 to 13 show the dose 
profiles for incident neon beams of energies 670 MeV / 
nucleon, 2 GeV/nucleon, and 10 GeV/nucleon; part a 
of each figure presents the individual ion doses and 
part b the normalized total dose D(x)/D( 0). In each 
case, the Bragg profile for the incident ion is clearly 
seen. As the incident energy increases, the maximum 
flux for the lighter ions moves farther into the shield 
and away from the Bragg cutoff at R Po /vj. This 
movement causes the total dose to extend past the 
Bragg cutoff. 

Unlike in the previous cases, when 5-GeV/nucleon 
cobalt ions are the incident ions (figs. 14(a) and 
14(b)), a significant amount of the dose occurs after 
the Bragg cutoff («105 g/cm 2 ). This is a result of the 
larger number of fragments produced and the deeper 
penetration of these fragments at a given energy per 
nucleon. Also, the larger number of ions gives rise 
to an increase in the dose near x = 0, rather than a 
decrease as seen for incident neon. 

The following dose calculations were performed 
with the more realistic stopping power given by 
equations (27) and (29a). The integrals in equa- 
tion (26c) were performed numerically by using a 
reliable Romberg integration routine (ref. 6). For in- 
cident neon of 670 MeV/nucleon, figure 15(a) shows 
the familiar Bragg peak and cutoff, and the contribu- 
tion of the lighter ions is shown in figure 15(b). The 
profile of ion j — 9 seems to follow the dose profile of 
the incident ion near the Bragg peak. This aspect of 
course results from the simple fact that the incident 
ion is the sole producer of ion j = 9. As j decreases, 


the maximum dose for each value of j moves far- 
ther into the medium away from the Bragg peak as 
the production of each subsequent ion becomes less 
dependent on the source ion. The similarity of the 
dose profiles for ions j = 9 and 10 is even more strik- 
ing for an incident energy of 1 GeV/nucleon. (See 
figs. 16(a) and (b).) For higher incident energies 
(figs. 17(a) and (b)), the Bragg peak dissipates as 
a result of the increased fragmentation at larger pen- 
etration distances. Figure 18 presents the total dose ■ 

for the previous three cases. The Bragg peaks again 
stand out for the lower incident energies. Figures 19 
and 20 show the doses for nickel incident at energies 
of 1 GeV/nucleon and 5 GeV/nucleon, respectively. j 

The curves represent only ions 18 < j < 28, because, I 

for j < 18, catastrophic round-off error accumulation I 

occurs and limits the evaluation. j 

i 

i 

As an example of the use of the analytical solu- j 

tion as a predictive tool, the experimental dose from a j 

neon beam incident on water (ref. 11) was compared ! 

with that obtained from the analytical solution. The | 

comparison of the normalized doses is shown in fig- j 

ure 21. The curves are normalized to the dose at 
x = 0, and cross-section normalization a was ad- f 

justed slightly to obtain a reasonable fit to the data. I 

If desired, the fragmentation parameter could also be 
adjusted to obtain an even better fit. [ 

[ 

4.3. Distributed Sources j 

£ 

Numerical results are now presented for analytical [ 

and solar distributed sources. An analytic source 1 

is one that allows the integral value of J r to be r 

expressed in closed form; the solar source requires 
a numerical quadrature evaluation. i 

4-3.1. Analytic source. Figures 22 and 23 
show the flux profiles for neon ions that are inci- i 

dent with a distributed source in s (or E ), given { 

by equation (37a) or 37(b), and a source cutoff 
taken at 10 4 MeV/nucleon. For small values of M 0 
(0.01) and an energy very near the cutoff (E = 

9999.9 MeV/nucleon), the incident ion source dom- 
inates all other sources from fragmentation. Thus, 
in the region where the incident ion has a nonzero 
flux (0 < x < sjvj ), all lighter ion flux profiles are 
essentially identical. For this energy, the ions have 
only penetrated ^0.01 g/cm 2 , and significant frag- 
mentation has not occurred. As the energy decreases 
(s increases) to 1 GeV/nucleon, the lighter ions build 
up similar to the beam source. 

Similar behavior is observed in the flux profiles 
for a source that is more peaked at s = 0 (M 0 ~ 

1.0). As shown in figure 23(a), the profiles are 
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virtually identical to the preceding case, since over 
the very short penetrations both sources look alike. 
For an energy of 1 GeV/nucleon (fig. 23b) the source 
resembles a delta function because of its rapid decay 
in s, and the profiles look like those for a beam source. 


(^200 g/cm 2 ), the neutron flux decays exponentially 
(not shown), since there is no longer a neutron source 
in this region — only absorption. 

4.6. Comparisons With Numerical 
Solutions 


4*3.2. Solar source . A source of the form 
f(E) — Qo/(b\ + E) n 

is typically taken to represent solar ions from the 
Sun. The J r integral in equation (37c) cannot be 
performed analytically, so the numerical Romberg 
integration must be used. (See ref. 12.) Because 
of the excessive computer time required for these 
calculations, only a few ions (<10) can be considered. 
Figures 24(a) and (b) show the resulting ion fluxes 
(to j = 19) as a function of E for incident nickel with 
Qo = 28 446, b\ = 1000, and n = 2.5 at positions 
x = 1 g/cm 2 and 10 g/cm 2 . As the distance into the 
shield increases, the lighter ion fluxes again increase 
because of fragmentation. All ion fluxes seem to have 
an energy variation similar to the incident ion source; 
this similarity suggests that a scaling law may apply. 

4.4. Composite Beam 

Figures 25(a) to 25(c) show the flux profiles using 
the Wilson stopping power for an artificial cosmic ray 
beam composed of the elements neon through iron 
with the assumed abundances given in table V. The 
jagged appearance of the flux profiles near the source 
energy 4.9 GeV/nucleon is a result of the relatively 
short penetration distances and the associated lack of 
fragmentation so close to the source energy. For this 
reason, the profiles for the ions with a charge near 
that of the incident ion still reflect the source distri- 
butions. The lighter ions, however, have smoother 
flux profiles as more fragmentation occurs at lower 
energies. 

4.5. Generation of Neutrons 

The parameter f n in the cross-section model 
(eq. (45)) has been introduced to make the distinc- 
tion between the neutron and proton components. 
(See section 2.6.) In figures 26(a) and (b), the re- 
sulting ion fluxes at E = 500 MeV/nucleon for an in- 
cident neon beam of E 0 = 1 GeV/nucleon are shown 
for fn = 0 and 0.5. In addition, the Wilson stop- 
ping power has been used. As expected when f n is 
increased from zero to one, the neutron component 
also increases with an attendant decrease in the pro- 
ton component. The neutron component does not 
go to zero at a finite distance as for charged par- 
ticles. After a distance equal to the proton range 


4-6.1 . Finite- difference /power series ex- 
pansion solution. For G intervals, if the GIT 
equations are discretized in the path-length variable, 
equations (3a) and (3b) become 


[b + (°i + £)] ^ (x) = £ *r ,(i) 


j 

+ H 03 ]k 4^ 

k—j+1 

) = fj 6jj (9= 1 , 2 , 


(47a) 
G) (47b) 


The advantage of the discretized model is that 
energy-dependent cross sections and multiplicities 
can now readily be included. The solution to equa- 
tions (47a) and (47b) is most easily obtained by ap- 
plying the power series method, where the expansion 

oo ^9^ 

4>j,g ( x ) = E, “T “ ar -^ ( 48a ) 

71=0 n ' 


is introduced into equations (47) to yield the follow- 
ing recursion relation for h g n T -\ 


hl r . = - L . + 2L) h «- T + 2l h 9-y ) 

n,j y Ji9 J n,j As n — l,j I 


+ E *r.« 

A:=j + 1 


(48b) 


The particular form of equation (48a) is valid for 
x in the interval [a r _i,a r ]. For incident neon at 
E 0 = 670 MeV/neutron, a flux comparison of the 
above interval solution with the analytical solution 
(eq. (22a)) is given in figure 27. The flux energy 
spectrum is shown at a position of 10 g/cm 2 from the 
surface. Even for relatively large values of As, the 
proton flux (j = 1) is in relatively good agreement; 
there is poor agreement for the heavy ions. This 
poor agreement is probably a result of the highly 
singular nature of the incident beam (eq. (4a)), which 
leads to discontinuous distributions for ions 9 and 10; 
these ions are poorly simulated with finite differences. 
However, as As is reduced, the general agreement 
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improves. Another area of disagreement is at the 
beginning and end of the ion travel, where the flux 
falls sharply to zero. Because of numerical diffusion, 
the finite-difference solution does not exhibit step 
function behavior, but the fluxes do become steeper 
with increasing As. The finite-difference solution 
does, however, give the correct flux average at these 
discontinuities as As decreases. 

4.6.2. Integral transport finite- difference 
solution . Our final comparison is with the GCR 
code (ref. 11) currently in use at Langley Research 
Center and Johnson Space Center. The GCR code 
solves the GIT equations by integration along the ion 
trajectories on a spatial and energy grid. From pre- 
vious benchmark comparisons, the code is expected 
to be accurate to better than 10 percent. 

Figure 28 shows the flux profiles for the GCR code 
and benchmark for a distributed source of neon of the 
form 

f{s) = exp(O.Ols)0(lO 5 - E) 

impinging on an air shield. The nuclear fragmen- 
tation model specified by equations (46) was used 
along with the simplified stopping power given by 
equations (29). To provide a contrast, the bench- 
mark calculation begins at E = 20 MeV, whereas 
the GCR calculation begins at E = 100 MeV. The 
two evaluations are in excellent agreement even at 
100 g/cm 2 . Figure 29 shows the flux profiles for 
incident nickel with the same source; these profiles 
again indicate virtually complete agreement. These 
comparisons provide added confidence that the GCR 
code employs a very accurate transport algorithm. 

5. Summary and Recommendations 

Based on its success in generating previous bench- 
mark solutions, a straightforward analytical ap- 
proach has been adopted to solve the full GIT 
equations. As a result, a closed-form analytical 
representation of the ion fluxes that result from a 
monocnergetic incident ion source has been derived. 
While there are no truncation errors to contend with, 
round-off error docs present a problem. The prob- 
lem becomes serious only when the incident ion is 
very heavy (J > 28) and when the energies are near 
the source energy; otherwise, 64-bit, double-precision 
arithmetic adequately allows the evaluation of the 


flux representation. The dose profiles can also be 
obtained from the analytical representation with the 
additional evaluation of an integral. 

Since the flux for the beam source is actually a 
Green’s function, the flux profiles for sources dis- 
tributed in energy can be obtained by integration 
over the source energy. This integration, unless it 
is performed analytically, introduces round-off error 
contamination that may limit the total number of 
isotope flux profiles that can be evaluated. A more 
general beam source composed of several ions with 
the same velocity can also be accommodated. 

Finally, a distinction between neutrons and pro- 
tons is made to more correctly take into account the 
attenuation of neutrons. The neutrons, however, are 
still monoenergetic, and therefore do not contribute 
to the overall dose. It is this last point that deserves 
further attention. 

An improved neutron flux model, including scat- 
tering interactions, should be the first step of a fu- 
ture effort. By using the ions where j > 3 to specify 
a neutron source, it is possible to formulate a more 
realistic Boltzmann equation for the neutron compo- 
nent. A multigroup treatment can then be applied 
to this equation both with and without the straight- 
ahead assumption. In this way, the neutron energy 
spectrum can be obtained either by a discretization 
in the spatial variable or a discrete-ordinates formu- 
lation in the angular variable. Once angular depen- 
dence has been introduced into the formulation, two- 
and three-dimensional shields can be considered. 

In addition to continuing to develop benchmarks 
for more realistic problems, such as including en- 
ergy dependence in the cross sections or angular de- 
pendence, the benchmark solutions developed thus 
far should be compared with existing transport al- 
gorithms to establish a common basis of compari- 
son for these algorithms. The benchmark solutions 
should be considered as a package that includes the 
energy-independent, spatially independent and cou- 
pled energy-space solutions. Each solution should be 
checked to verify each specific part of a transport 
algorithm. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
August 26, 1991 
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7. Symbols 

r 0 

defined in equation (26d) 

Aj 

mass number of jth ion 

S 3 

stopping power of jth ion, 

dj 

position of zth interval boundary, 
g/cm 2 

Sv 

MeV-cm 2 /g 

proton stopping power, 

a £m 

a\ 

b(m 

Of 

= &£ ~ &m 

= oj - 1 - 
= V£- Vm 

s(E, Eg) 

X 

MeV-cm 2 /g 

proton path length, g/cm 2 
position, g/cm 2 

h 

III 

£ 

1 

i 

y'k 

ij 

defined by equation (7b) 

Dj 

dose due to jth ion, MeV/g 

ij 

defined by equation (22b) 

dj 

abundance of jth beam ion 

a 

a ij 

Ol 0) Qq , Tl 0 

_ l 

E 

fj 

energy per nucleon, MeV/nucleon 
monoenergetic beam source term 

_ Ot\OL 0 

recursion coefficient 
parameters used in equation (27) 

f 3 
h 9 ’ r 

n, ] 

Ir 

Laplace transform of jj 
power series coefficient 

(3 

_ _i_ 
— oc 0 

defined by equation (26c) 

Sjk 

= M jk°k (eq- (3c)) 

7r,i 

defined by equation (31a) 

As 

finite-difference interval, g/cm 2 

J 

charge number of incident ion 

6(s) 

Dirac 6-function 

3 f 

defined by equation (36b) 

&iJ 

Kronecker delta 

j 

charge number of jth ion 

6(x) 

Heaviside step function 

h 

defined by equation (40a) 

Af 

defined by equation on page 4 

Lf 

defined by equation (39b) 

U 3 

= jfc ( e< 4- (2b)) 

Mjk 

multiplicity of ion j produced by 
ion k 

a 

j 

cross-section parameter, cm 2 /g 

N e 

defined by equation (12) 

a 3 

macroscopic absorption cross 

N p 

number of intervals 


section for jth ion, cm 2 /g 

Q(x,s) 

volume source coefficient 

<fj 

flux of jth ion 

Qo ) Mo 

parameters in equation (37a) 

f)j 

Laplace transform of 4>j 

R 

range, g/cm 2 

X 

defined by equation on page 5 
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to t> oo cn> 


Table I. Relation of r Index to £ 

A r (p) — a mn + bmnp\ 


e 

m 

n 

N e 

1 

j - 1 

j 

1 

2 

J — 2 

j 

3 


J- 2 

j - i 


3 

J - 3 

j 

6 


J - 3 

j - i 



J — 3 

J — 2 


4 

J — 4 

J 

10 


J — 4 

J- 1 



J — 4 

J — 2 



J — 4 

J- 3 



Table II. Variation of i\ and Z 2 To Satisfy Equation (20) With r 


i 

T 

n 

*2 

1 

1 

0 

i 

2 

2 ' 

0 

2 


3 

1 

2 

3 

4 

0 

3 


5 

1 

3 


6 

2 

3 

4 

7 

0 

4 


8 

1 

4 


9 

2 

'4 


10 

3 

4 

5 

11 

0 

5 


12 

1 

5 


13 

2 

5 


14 

3 

5 


15 

4 

5 






Tabic IV. Partial Sums in Evaluation of Equation (22a) 

m 

S m = E T r 

r= 1 


r 

T r 

Sm 

r 

T r 

Sm 

37 

3.9228E+07 

3.9228E+07 

79 

— 1.1811E+08 

— 7.5948E+10 

38 

9.7776E+08 

1.0170E+09 

80 

— 8.1489E+09 

-8.4097E+10 

39 

1.1020E+10 

1.2037E+10 

81 

1.2471E+10 

— 7.1626E+10 

40 

— 1.5280E+06 

1.2035E+10 

82 

1.6089E+11 

8.9259E+10 

41 

8.7956E+06 

1.2044E+10 

83 

1.6572E+10 

1.0583E+11 

42 

-1.5549E+02 

1.2044E+10 

84 

3.8245E+05 

1.0583E+11 

43 

1.3975E+02 

1.2044E+10 

85 

2.2172E+08 

1.0605E+11 

44 

1.2963E+01 

1.2044E+10 

86 

— 1.1062E+11 

-4.5652E+09 

45 

1.7750E+01 

1.2044E+10 

87 

1.2098E+06 

— 4.5640E+09 

46 

— 1.7809E+03 

1.2044E+10 

92 

2.2978E+10 

1.8414E+10 

47 

— 3.3768E+04 

1.2044E+10 

93 

— 1.8297E+09 

1.6584E+10 

48 

1.6875E+04 

1.2044E+10 

94 

— 7.1510E+09 

9.4334E+09 

49 

2.2094E+03 

1.2044E+10 

95 

— 1.6057E+08 

9.2729E+09 

50 

— 1.1694E+05 

1.2044E+10 

96 

3.5538E+08 

9.6282E+09 

51 

— 5.3195E+01 

1.2044E+10 

97 

8.8145E+04 

9.6283E+09 

52 

1.3282E+05 

1.2044E+10 

98 

— 1.7218E+06 

9.6266E+09 

53 

6.4344E+02 

1.2044E+10 

99 

2.4350E+07 

9.6510E+09 

54 

— 3.7529E— 02 

1.2044E+10 

100 

— 1.4887E+08 

9.5021E+09 

56 

— 3.8850E+06 

1.2040E+10 

106 

1.3371E+08 

9.6358E+09 

57 

— 5.6258E+09 

6.4145E+09 

107 

1.2474E+10 

2.2110E+10 

58 

2.7668E+08 

6.6911E+09 

108 

— 1.1917E+10 

1.0192E+10 

59 

-1.1491E+11 

— 1.0821E+11 

109 

8.1157E+08 

1.1004E+10 

60 

— 2.0332E+10 

— 1.2855E+11 

110 

5.9305E+09 

1.6934E+10 

61 

— 1.7720E+05 

— 1.2855E+11 

111 

2.4703E+04 

1.6934E+10 

62 

1.9891E+07 

— 1.2853E+11 

112 

— 4.6283E+06 

1.6930E+10 

63 

6.6318E+10 

— 6.2209E+10 

113 

— 2.2498E+09 

1.4680E+10 

64 

— 5.9047E+03 

— 6.2209E+10 

114 

— 1.3697E+07 

1.4666E+10 

67 

-1.3597E+10 

-7.5806E+10 

121 

— 9.4346E+09 

5.2317E+09 

68 

— 2.5725E+08 

— 7.6063E+10 

122 

2.3988E+09 

7.6304E+09 

69 

1.8921E+08 

— 7.5874E+10 

123 

— 4.8884ET09 

2.7420E+09 

70 

1.0403E+07 

-7.5864E+10 

124 

6.2213E+09 

8.9633E+09 

71 

— 3.9645E+07 

— 7.5903E+10 

125 

— 2.4306E+09 

6.5326E+09 

72 

— 5.2003E+02 

-7.5903E+10 

126 

— 3.1752E+05 

6.5323E+09 

73 

1.1232E+04 

— 7.5903E+10 

127 

— 2.3539E+08 

6.2969E+09 

74 

— 7.1222E+04 

— 7.5903E+10 

128 

— 6.3879E+09 

— 9.0967E+07 

75 

7.3132E+07 

-7.5830E+10 

129 

9.0967E+07 

9.7023E— 04 
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Table V. Composite Cosmic Ray Beam 


Element 

Abundance, percent 

Fe 

1.3 

Mn 

0.6 

Cr 

1.7 

V 

1.4 

Ti 

2.0 

Sc 

2.0 

Ca 

1.8 

K 

2.5 

Ar 

1.8 

Cl 

2.6 

s 

2.5 

p 

5.1 

Si 

7.2 

A1 

11.5 

Mg 

16.3 

Na 

14.7 

Ne 

23.1 

Fet 

1.9 


* Portion of cosmic ray spectrum from neon 
to iron. 


t Added to make 100 percent. 
















(a) a = 10 -4 cm 2 /g (rarefied shield). 



(b) a = 1 cm 2 /g (dense shield). 

Figure 9. Sensitivity of flux at 995 MeV/nucleon to incident neon (J = 10) ions at E 0 = 1 GeV/nucleon with 
cross-sect io n normalization . 







(b) Sum of all ions. 


Figure 12. Dose profiles ( n 0 = 1) for incident neon (J — 10) ions at E 0 — 2 GeV/nucleon. 





(b) Sum of all ions. 

Figure 14. Dose profiles (n 0 — 1) for incident cobalt (J ~ 27) ions at E 0 = 5 GeV. 




>0 



(a) J = 10. 



(b) j < 10. 

Figure 15. Dose profiles ( n 0 = 1.79) for incident neon (J = 10) ions at E () = 670 MeV/nucleon 
















Figure 20. Dose profiles ( n 0 = 1.79) for incident nickel (J = 28) ions at E 0 — 5 GeV/nucleon. 



Figure 21. 


Experimental dose comparison for incident neon ions 


(J = 10) on water (cr ^ 1.5 x 10 2 


cm 2 /g) 
















(a) E = 4.9 GeV/nucleon. 



(b) E = 4.3 GeV/nucleon. 

Figure 25. Flux profiles for a composite beam of neon through iron ions at E 0 = 5 GeV/nucleon. 
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50 75 100 125 150 175 200 225 

x, g/cm 2 

(a) /„ = 0 (no neutrons). 



(b) fn = 0.5. 

Figure 26. Distinction between neutrons and protons at E = 500 MeV/nucleon for incident neon ( J = 10) ions 
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